load Figure4s_matlab_rawdata.mat
%% Map to real gate voltage
Vg=Vg1;

pmin = -6;
pmax = -2;
nmin = 5;
nmax = 9;
vg=Vg;
vg(vg<0.5) = pmin + (pmax-pmin)*vg(vg<0.5)*2;
vg(vg>=0.5) = nmin + (nmax-nmin)*(vg(vg>=0.5)*2-1);

% Magnet axis
th0 = 106/180*pi;
%%

figure(4967);
clf;

figure(4968);
clf;

tm = zeros(size(Rxx1,1),6);
amp = zeros(size(Rxx1,1),6);

for k=1:6
    switch(k)
        case 1
            Rxx=Rxx1;
            Bpar=squeeze(Bpar1(1,1,:));
            Btheta=Btheta1;
        case 2
            Rxx=Rxx2;
            Bpar=squeeze(Bpar2(1,1,:));
            Btheta=Btheta2;
        case 3
            Rxx=Rxx3;
            Bpar=squeeze(Bpar3(1,1,:));
            Btheta=Btheta3;
        case 4
            Rxx=Rxx4;
            Bpar=squeeze(Bpar4(1,1,:));
            Btheta=Btheta4;
        case 5
            Rxx=Rxx5;
            Bpar=squeeze(Bpar5(1,1,:));
            Btheta=Btheta5;
        case 6
            Rxx=Rxx6;
            Bpar=squeeze(Bpar6(1,1,:));
            Btheta=Btheta6;
    end

    Rxx(isnan(Rxx))=0;
    Rxx(isinf(Rxx))=0;
    TH=squeeze(Btheta(1,:,1));
    G=vg(:,1);

    [~,ind] = min(abs(Bpar-1));
    R = squeeze(Rxx(:,:,ind));
    A = R ./ repmat(max(R,[],2), 1, size(R,2));
    
    figure(4967);
    subplot(2,3,k);
    surf(TH-th0, G, A, 'EdgeColor','none'); view(2); axis tight;
    caxis([0,1]);
    ylabel 'Vg'
    xlabel '\theta'
    
    figure(4968);
        
    th = linspace(0,2*pi,size(A,2))';
    s = sin(2*th);
    c = cos(2*th);
    
    for i=1:size(A,1)
        a=A(i,:);
        f = fft(a);
        % Get phase of the 2theta component
        tm(i,k) = mod(pi/2-angle(f(3))/2-th0, pi);
        amp(i,k)= (max(a)-min(a)) / (max(a)+min(a));
        if amp(i,k)<0 || amp(i,k)>1
            amp(i,k)=NaN;
            tm(i,k)=NaN;
        end
       
    end
    
    subplot(1,2,1);
    hold on
    plot(G,(mod(tm(:,k)+pi/2,pi)-pi/2)/pi*180, '-o', 'Color', [k*0.16,0.2,0.1]);
    ylim([-90,90]);
    ylabel 'Anisotropy angle (deg)';
    xlabel 'V_g';
    
    subplot(1,2,2);
    hold on
    plot(G, amp(:,k),'-o', 'Color', [k*0.16,0.2,0.1]);
    ylim([0,1]);
    ylabel 'Anisotropy amplitude';
    
end

%% Export data

fp = fopen('DeviceB-anisotropy-ampphase.txt', 'w');
fprintf(fp, '# Vg A(base) th(base) A(0.4K) th(0.4K) A(0.9K) th(0.9K) A(1.7K) th(1.7K) A(2.4K) th(2.4K) A(4K) th(4K)\n');
fprintf(fp,[repmat('%g\t', 1, 13), '\n'], [vg(:,1,1), amp(:,1), tm(:,1), amp(:,2), tm(:,2), amp(:,3), tm(:,3), amp(:,4), tm(:,4), amp(:,5), tm(:,5), amp(:,6), tm(:,6)]');
fclose(fp);
